Monomer-dimer model in two-dimensional rectangular lattices with fixed dimer 

density 
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The classical monomer-dimer model in two-dimensional lattices has been shown to belong to 
the "#P-complete" class, which indicates the problem is computationally "intractable". We use 
exact computational method to investigate the number of ways to arrange dimers on m x n two- 
dimensional rectangular lattice strips with fixed dimer density p. For any dimer density < p < 1 , 
we find a logarithmic correction term in the finite-size correction of the free energy per lattice site. 
The coefficient of the logarithmic correction term is exactly —1/2. This logarithmic correction term 
is explained by the newly developed asymptotic theory of Pemantle and Wilson. The sequence 
of the free energy of lattice strips with cylinder boundary condition converges so fast that very 
accurate free energy /2(p) for large lattices can be obtained. For example, for a half-filled lattice, 
/ 2 (l/2) = 0.633195588930, while / 2 (l/4) = 0.4413453753046 and / 2 (3/4) = 0.64039026. For p < 
0.65, / 2 (p) is accurate at least to 10 decimal digits. The function /2(p) reaches the maximum value 
/2(p*) = 0.662798972834 at p* = 0.6381231, with 11 correct digits. This is also the monomer-dimer 
constant for two-dimensional rectangular lattices. The asymptotic expressions of free energy near 
close packing are investigated for finite and infinite lattice widths. For lattices with finite width, 
dependence on the parity of the lattice width is found. For infinite lattices, the data support the 
functional form obtained previously through series expansions. 

PACS numbers: 05.50.+q, 02.10.Ox, 02.70.-c 



I. INTRODUCTION 

The monomer-dimer problem has received much atten- 
tion not only from statistical physics but also from the- 
oretical computer science. As one of the classical lattice 
statistical mechanics models, the monomer-dimer model 
was first used to describe the absorption of a binary mix- 
ture of molecules of unequal sizes on crystal surface . 
In the model, the regular lattice sites are either covered 
by monomers or dimers . The diatomic molecules are 
modeled as rigid dimers which occupy two adjacent sites 
in a regular lattice and no lattice site is covered by more 
than one dimer. The lattice sites that are not covered 
by the dimers are regarded as occupied by monomers. A 
central problem of the model is to enumerate the dimer 
configurations on the lattice. In 1961 an elegant ana- 
lytical solution was found for a special case of the prob- 
lem, namely when the planar lattice is completely covered 
by dimers (the close-packed dimer problem, or dimer- 
covering problem) 0, Q ■ For the general monomer-dimer 
problem where there are vacancies (monomers) in the 
lattice, there is no exact solution. For three-dimensional 
lattices, there is even no exact solution for the special 
case of close-packed dimer problem. One recent advance 
is an analytic solution to the special case of the prob- 
lem in two-dimensional lattices where there is a single 
vacancy at certain specific sites on the boundary of the 
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lattice 0, Q . The monomer-dimer problem also serves 
as a prototypical problem in the field of computational 
complexity 6J. It has been shown that two-dimensional 
monomer-dimer problem belongs to the '#P- complete" 
class and hence is computationally intractable 0. 

Even though there is a lack of progress in the analytical 
solution to the monomer-dimer problem, many rigorous 
results exist, such as series expansions M, M, Hfj| . lower 
and upper bounds on free energy [Til Il2|. monomer- 
monomer correlation function of two monomers in a lat- 
tice otherwise packed with dimers 01 > locations of zeros 
of partition functions, [13. Il5j. and finite-size correction 
[l6j . Some approximate methods have also been pro- 
posed [l7L llcl Il9l l2Cj . The monomer-dimer constant 
hd (the exponential growth rate) of the number of all 
configurations with different number of dimers has also 
been calculated 0, 0] . By using sequential importance 
sampling Monte Carlo method, the dimer covering con- 
stant for a three-dimensional cubic lattice has been es- 
timated |2l| . The importance of monomer-dimer model 
also comes from the fact that there is one to one map- 
ping between the Ising model and the monomer-dimer 
model: the Ising model in the absence of an external 
field is mapped to the pure dimer model [22, [U 0, HH j 
and the Ising model in the presence of an external field 
is mapped to the general monomer-dimer model [l4| . 

The major purposes of this paper are (1) to show it 
is possible to calculate accurately the free energy of the 
monomer-dimer problem in two-dimensional rectangular 
lattices at a fixed dimer density by using the proposed 
computational methods (Sections ITT1 IIVI IVII and IVIl(l . 
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and (2) to use the computational methods to probe the 
physical properties of the monomer-dimer model, espe- 
cially at the high dimer density limit ( Section I VIII jl . The 
high dimer density limit is considered to be more dif- 
ficult and more interesting than the low dimer density 
limit. The major result is the asymptotic expression Eq. 
1241 The third purpose of the paper is to introduce the 
asymptotic theory of Pemantle and Wilson [26j, which 
not only gives a theoretical explanation of the origin of 
the logarithmic correction term found by computational 
methods reported in this paper fSection llllfl . but also has 
the potential to be applicable to other statistical models. 

The following notation and definitions will be used 
throughout the paper. The configurational grand canon- 
ical partition function of the monomer-dimer system in 
a m x n two-dimensional lattice is 



^mn(i) =ajv(m,fi)i +aAr_i(m, n)x +■ • -+ao(m, n) 

(1) 

where a s (m, n) is the number of distinct ways to arrange 
s dimers on the m x n lattice, N — [mn/2\, and x can 
be taken as the activity of a dimer. The average number 
of sites covered by dimers (twice the average number of 
dimers) of this grand canonical ensemble is given by 



2 d\nZ m ^ n (x) _ 2 2 s =i a s{m, n)sx 



d In: 



EfLo a s (m,n)x s 



(2) 

The limit of this average for large lattices is denoted as 
9(x): 9(x) — lim m , n ->oo 8 m ,n{x). In general we use Od(x) 
for the average number of sites covered by dimers in a 
d-dimcnsional infinite lattice when the dimer activity is 
x. 

The total number of configurations of dimers is given 
by Z m>n (V) at x = 1, and the monomer-dimer constant 
for a two-dimensional infinite lattice is defined as 



h-r = lim 



hiZ m , n (l) 



m.n — >oo 



(3) 



In general, we denote hd as the monomer-dimer constant 
for a ci-dimensional infinite lattice, and hd(x) as the grand 
potential per lattice site at any dimer activity x. For a 
two-dimensional infinite lattice, 



hzix) = lim 



m.n — >oo 



(4) 



In this paper we focus on the number of dimer con- 
figurations at a given dimer density p. In this sense we 
are working on the canonical ensemble. The connection 
between the canonical ensemble and the grand canoni- 
cal ensemble is discussed in Appendix [S] We define the 
dimer density for the canonical ensemble as the ratio 



P : 



2s 
mn 



(5) 



When the lattice is fully covered by dimers, p = 1. For 
a m x n lattice, the number of dimers at a given dimer 



density is s 



mnp 



In the following we use a m ,n(p) as the 



number of distinct dimer and monomer configurations at 
the given dimer density p. By using this definition, Eq. 
^can be rewritten as 



Zm,n(%) — ^ ' ®m,nip)^ ^ ■ 



(6) 



0<p<l 



The free energy per lattice site at a given dimer density 
p is defined as 



fm,n(p) 



In a m ,n(p) 



and the free energy at a given dimer density for a semi- 
infinite lattice strip oo x n is 

/oo.n(p) = lim ln n "'.»(p) = lim f ( p y 

m^oo 77777 m— >oo 

For infinite lattices where both m and n go to infinity, 
the free energy is 

h{p) = /oo.co(p) = lim ln = lim f Q,y 

m,n — >oo 77777 n — >oo 

We use the subscript d in fd(p) to indicate the dimension 
of the infinitely large lattice. From the exact result @>01 
we know h(p) at p = 1 



Mi) 



G 



0.291560904 



where G is the Catalan's constant. For other values of 
p > 0, no analytical result is known, although several 
bounds are developed [Hill- 

We will show below that 
by using the exact calculation method developed previ- 
ously |23, 123, 123 1 we can calculate / 2 (ft) at an arbitrary 
dimer density p with high accuracy. 

The article is organized as follows. In Section [HI the 
computational method is outlined. In Section IIIII we 
show a logarithmic correction term in the finite-size cor- 
rection of f m ,n(p) for any fixed dimer density < p < 1. 
The coefficient of this logarithmic correction term is ex- 
actly — 1/2, for both cylinder lattices and lattices with 
free boundaries. We give a theoretical explanation for 
this logarithmic correction term and its coefficient us- 
ing the newl y d eveloped asymptotic theory of Pemantle 
and Wilson |2|j. In this section we point out the uni- 
versality of this logarithmic correction term with coeffi- 
cient of —1/2. This term is not unique to the monomer- 
dimer model: a large class of lattice models has this 
term when the "density" of the models is fixed. More 
discussions of applications of this asymptotic method to 
the monomer-dimer model in particular, and statistical 
models in general, can be found in Section IIXI In Sec- 
tion ]W\ we calculate foo,n(p) on lattice strips 00 x n for 
77 = 1, . . . , 17 with cylinder boundary condition. The 
sequence of /oo,n(p) on cylinder lattices converges very 
fast so that we can obtain h{p) quite accurately. To 
the best of our knowledge, the results presented here are 
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the most accurate for monomer-dimer problem in two- 
dimensional rectangular lattices at an arbitrary dimer 
density. In Section similar calculations of foc,n{p) are 
carried out on lattice strips oo x n with free boundaries 
for n = 1 , . . . , 16. Compared with the sequence with 
cylinder boundary condition, the sequence foa,n(p) with 
free boundaries converges slower. In Section IVII the po- 
sition and values of the maximum of f2{p) are located: 
f 2 (p*) = 0.662798972834 at p* = 0.6381231. These re- 
sults give an estimation of the monomer-dimer constant 
with 11 correct digits. The previous best result is with 
9 correct digits [lj- The results are also compared with 
those obtained by series expansions and field theoreti- 
cal methods. The maximum value of f2(p) is equal to 
the two-dimensional monomer-dimer constant h,2- This 
is one special case of the more general relations between 
the calculated values in the canonical ensemble and those 
in the grand canonical ensemble, and these relations are 
further discussed in Appendix 1X1 In Section lYlII the re- 
lations developed in Appendix[X]are used to compare the 
results of the computational method presented in this pa- 
per with those of Baxter jl9j | . For monomer-dimer model, 
the more interesting properties are at the more difficult 
high dimer density limit. In Section [Villi asymptotic be- 
havior of the free energy foo,n(p) is examined for high 
dimer density near close packing. For lattices with finite 
width, a dependence of the free energy foo,n(p) on the 
parity of the lattice width n is found (Eq. 1221) , consistent 
with the previous results when the number of monomers 
is fixed [29j . The combination of the results in this sec- 
tion and those of Section ITTT1 leads to the asymptotic ex- 
pression Eq. [2] for near close packing dimer density. 
The asymptotic expression of f2(p), the free energy on 
an infinite lattice, is also investigated near close packing. 
The results support the functional forms obtained pre- 
viously through series expansions 0, but quantitatively 
the value of the exponent is lower than previously con- 
jectured. In Appendix [E] we put together in one place 
various explicit formulas for the one-dimensional lattices 
(n = 1). These formulas can be used to check the for- 
mulas developed for the more general situations where 
n > 1. As an illustration, an explicit application of the 
Pemantle and Wilson asymptotic method is also given 
for n = 1. 



II. COMPUTATIONAL METHODS 

The basic computational strategy is to use exact calcu- 
lations to obtain a series of partition functions Z m ^ n (x) 
of lattice strips m x n. Then for a given dimer density p, 
fm,n{p) can be calculated using arbitrary precision arith- 
metic. By fitting f m ,n(p) to a given function (Sections 
HVllVl and lVIIl|) . /oo,n(p) can be estimated with high ac- 
curacy. From /oo,n(/o), fa(,p) can then be estimated using 
the special convergent properties of the sequence foo,n{p) 
on the cylinder lattice strips ( Sect ion llV|) . 



A. Calculation of the partition functions 

The computational methods used here have been de- 
scribed in details previously [13, IH, ^ ■ The full par- 
tition functions (Eq. ^) are calculated recursively for 
lattice strips on cylinder lattices and lattices with free 
boundaries. As before, all calculations of the terms 
a s (m, n) in the partition functions use exact integers, and 
when logarithm is taken to calculate free energy fm.n(p), 
the calculations are done with precisions much higher 
than the machine floating-point precision. The bignum 
library used is GNU MP library (GMP) for arbitrary 
precision arithmetic (version 4.2) [3fJ. The details of the 
calculations on lattices with free boundaries can be found 
in Ref. 28, so in the following only information on cylin- 
der lattices is given. 

For a m x n lattice strip, a square matrix M n is set 
up based on two rows of the lattice strip with proper 
boundary conditions. The vector f2 m , which consists of 
partition function of Eq. ^ as well as other contracted 
partition functions [27| , is calculated by the following re- 
currence 

Q m = M n fi m _i. (7) 

Similar recursive method has also been used for other 
combinatorial problems, such as calculation of the num- 
ber of independent sets [31j. For a cylinder lattice strip, 
the matrix M n is constructed in a similar way as that 
with free boundaries |28(. The total valid number (v c (n)) 
and unique number (u c (n)) of configurations are given 
respectively by the generating function ^ n v c (n)x n — 
i(3 + x - x 3 )/(l - 3x - x 2 )/{x - l)/(x + 1) and the for- 
mula 

/ _ V"> ^{d)2 n / d IV"- 1 )/ 2 if n odd 

u c[n) - 2^ 2n + 1 2"/ 2 - 1 + 2"/ 2 - 2 if n even 

d\n \ 

where <p(m) is Euler's totient function, which gives the 
number of integers relatively prime to integer m. The 
size of matrix M n is u c (n) x u c (n). The first 17 terms 
of the sequence v c (n) are: 3, 10, 36, 118, 393, 1297, 
4287, 14158, 46764, 154450, 510117, 1684801, 5564523, 
18378370, 60699636, 200477278, and 662131473. The 
first 17 terms of the sequence u c {n) are: 2, 3, 4, 6, 8, 13, 
18, 30, 46, 78, 126, 224, 380, 687, 1224, 2250, and 4112. 
It is noted that the sequence u c (n) is exactly the same as 
that shown in column 2, Table 1 of Ref. ^2- Calculations 
based on the dominant eigenvalues of the matrices of the 
cylinder lattice strips for n — 4, 6, 8, and 10 are carried 
out by Runnels (3^1 ■ The sizes of M n for cylinder lattice 
strips are smaller when compared with the correspond- 
ing numbers for lattice strips with free boundaries |28| . 
which allows for calculations on wider lattice strips. For 
cylinder lattice strips, full partition functions are calcu- 
lated for n = 1, . . . , 17, with length up to m — 1000 for 
n = 1, . . . , 13, m = 880 for n = 14, m = 669 for n = 15, 
m = 474 for n = 16, and m = 325 for n — 17. The corre- 
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sponding numbers for lattice strips with free boundaries 
are reported in Ref. I281 



B. Interpolation for arbitrary dimer density p 

In this paper the main quantity to be calculated is 
/2(p). The starting point of the calculations is the full 
partition function Eq. ^for different values of n and to. 
Finite values of n and to only lead to discrete values of 
dimer density p, as defined in Eq. |SJ For example, when 
n = 7 and m = 11, the number of dimers s takes the val- 
ues of 0, 1, ... , 38, and dimer density p of this lattice can 
only be one of the following values: 0, 2/77, . . . , 76/77. In 
general, for fixed finite m and n, p can only be a rational 
number: 



P 

P = ~ 

q 

where p and q are positive integers. When p is expressed 
as a rational number, the number of dimers is given by 



mnp mnp 
2 = 2q ' 



This expression is only meaningful if mnp can be divided 
by 2q. When we write the grand canonical partition func- 
tion in the form of Eq. EJfor finite m and n, we implicitly 
imply that Eq. [5] is satisfied. 

In the following we use the rational dimer density 
p = p/q whenever possible so that the value of a m , n (p) 
can be read directly from the partition function of to x n 
lattices. Depending on the values of p and q, some dimer 
densities, such as p — 1/2, can be realized in many lat- 
tices, while others can only be realized in small number 
of lattices with special combinations of values of m and 
n. In many situations it becomes impossible to use ra- 
tional p. For example, in Section IvTl the location of the 
maximum of f2{p) is searched within a very small re- 
gion of p, and in Section IVIII in order to compare the 
results from different methods, p takes the output values 
of other computational methods • I n such situations, 
if the rational form of p were used, p and q would be- 
come so big that not enough data points which satisfy 
Eq. |S] could be found for the fitting in the to x n lattice 
strip. To calculate f m ,n(p) for an arbitrary real number 
p (0 < p < 1), interpolation of the exact data points is 
needed. Since full partition functions have been calcu- 
lated for fairly long lattice strips, proper interpolation 
procedure can yield highly accurate values of f m ,n(p) for 
an arbitrary real number p. For interpolation, we use 
the standard Bulirsch-Stoer rational function interpola- 
tion method (3^, 03 • F° r an Y rea l number p, Eq. [S] is 
used to calculate the corresponding number of dimers s, 
which may not be an integer. On each side of this value of 
s, 30 exact values of a s (m, n) are used (if possible) in the 
interpolation. If on one side there are not enough exact 
data points of a s (m,n), extra data points on the other 
side of s are used to make the total number of exact data 



points as 60. For the high dimer density case (Section 
I Villi . the total number of data points used is changed 
to 30. We also take care that no extrapolation is used: if 
p is greater than the maximum dimer density for a given 
mxn lattice, the data point from this lattice is not used. 
Let's look at the above example of 11 x 7 lattice again. 
For this lattice, the highest dimer density is 76/77. If 
calculation is done for a dimer density p = 0.99, since 
p = 0.99 is greater than 76/77 ~ 0.987, the data point 
from this lattice will not be used in the following steps 
to avoid inaccuracy introduced by unreliable extrapola- 
tions. 



C. Fitting procedure 

The fitting experiments are carried out by using the 
"fit" function of software GNUPLOT (version 4.0) 
on a 64-bit Linux system. The fit algorithm imple- 
mented is the nonlinear least-squares (NLLS) Levenberg- 
Marquardt method 36]. All fitting experiments use the 
default value 1 as the initial value for each parameter, and 
(8) each fitting experiment is done independently. As done 



previously |2ll29|. only those a m ^ n (p) with to > 100 are 
used in the fitting. Since a rn ^ n (p) is calculated for rela- 
tively long lattice strips (in the to direction, see Section 
111 A|) . the estimates of foo,n(p) are usually quite accu- 
rate, up to 12 or 13 decimal place. The accuracy for this 
fitting step is limited by the machine floating-point pre- 
cision, since GNUPLOT uses machine floating-point repre- 
sentations, instead of arbitrary precision arithmetic. We 
would have used the GMP library to implement a fitting 
program with arbitrary precision arithmetic. This would 
increase the accuracy in the estimation of /2(p) when 
p is small. For the major objective of this paper, i.e., 
to investigate the behavior of /2(p) when p — > 1 (Sec- 
tion [^nj, however, the current accuracy is adequate. 
At high dimer density limit, the convergence of foo.n(p) 
towards f2{p) is much slower than at low dimer density 
limit. With lattice width n < 17 used for the current cal- 
culations, foo,n{p) is far from converging to the machine 
floating-point precision when p — > 1. 
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III. LOGARITHMIC CORRECTIONS OF THE 
FREE ENERGY AT FIXED DIMER DENSITY 

For lattice strips mxn with a fixed width n and a given 
dimer density p, the coefficients a m . n {p) of the partition 
functions are extracted to fit the following function: 



In a m ^ n (p) 



mu m ,„^; Ci C2 
Jm,n{P) = = c tH 1 n 

mn m m m° to* n m 

(9) 

where c = /oo,n(p)- 

For both cylinder lattices and lattices with free bound- 
aries, the fitting experiments clearly show that I = —1/2, 
accurate up to at least six decimal place, for any dimer 



C3 C4 £ ln(m + 1) 

TO TO* 
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density < p < 1. This result holds for both odd n and 
even n. This is in contrast with the results reported ear- 
lier for the situation with a fixed number of monomers 
(or vacancies), where the logarithmic correction coeffi- 
cient depends on the number of monomers present and 
the parity of the width of the lattice strip [2Sl Efll ] . We 
notice that a coefficient —1/2 also appears in the loga- 
rithmic correction term of the free energy studied in Ref. 
U which is a special case of the monomer-dimer problem 
where there is a single vacancy at certain specific sites on 
the boundary of the lattice. 

For the general monomer-dimer model, to our best 
knowledge, this logarithmic correction term with coeffi- 
cient of exactly —1/2 has not been reported before in the 
literature. The recently developed multivariate asymp- 
totic theory by Pemantle and Wilson 26], however, gives 
an explanation of this term and its coefficient. This the- 
ory applies to combinatorial problems when the multi- 
variate generating function of the model is known. For 
univariate generating functions, asymptotic methods are 
well known and have been used for a long time. The 
situation is quite different for multivariate generating 
functions. Until recently, techniques to get asymptotic 
expressions from multivariate generating functions were 
"almost entirely missing" (for review, see Ref. 126"). The 



newly developed Pemantle and Wilson method applies to 
a large class of multivariate generating functions in a sys- 
tematic way. In general the theory applies to generating 
functions with multiple variables, and for the bivariate 
case that we are interested in here, the generating func- 
tion of two variables takes the form 



F(x,y) 



G(x,y) 
H(x,y) 



a sm x s y m (10) 

s.m— 



where G(x,y) and H(x,y) are analytic, and H(0, 0) ^ 
0. In this case, Pemantle and Wilson method gives the 
asymptotic expression as 



Gix^jjo) l -y H y (x ,y ) 

»sm ~ 7== Xq y Q t —. r [11) 

V27T y mQ{x a ,y a ) 

where (xq, yo) is the positive solution to the two equations 



. N dH dH 

H(x,y)=0, mx— = sy— (12) 



and Q(x,y) is defined as 



~(xH x )(yH y ) 2 - (yH y ){xH x ) 2 - [(yH y ) 2 (x 2 H xx ) + {xH x f {y 2 H yy ) - 2(xH x )(yH y )(xyH xy )}. 

I 



Here H x , H y , etc. are partial derivatives dH/dx, dH/dy, 
and so on. One of the advantages of the method over 
previous ones is that the convergence of Eq. Illl is uniform 
when s/m and m/s are bounded. 

For the monomer-dimer model discussed here, with n 
as the finite width of the lattice strip, m as the length, 
and s as the number of dimers, we can construct the 
bivariate generating function F(x, y) as 

oo ron/2 oo 

F(x,y) =J2Y, as(m,n)x s y m = £ Z m . n (x)y m . 

7n—0 s—0 m— 

(13) 

For the monomer-dimer model, as well as a large class of 
lattice models in statistical physics, the bivariate gener- 



ating function F(x,y) is always in the form of Eq. i|10fl . 
with G(x,y) and H(x,y) as polynomials in x and y. In 
fact, we can get H(x,y) directly from matrix M n in Eq. 
(JZJ). It is closely related to the characteristic function 
of M n 27]: H(x,y) = det(M„ - I/y) x y w , where w is 
the size of the matrix M n . As an illustration, the bivari- 
ate generating function F(x, y) for the one-dimensional 
lattice (n = 1) is shown in Eq. IB 101 of Appendix iBl 

When the dimer density is fixed, which is the case dis- 
cussed here, s = pmn/2. If we substitute this relation 
into Eq. l(T2l . then we see that the solution (xo,yo) of 
Eq. (fT2"|) only depends on p and n, and does not depend 
on m or s. Substituting this solution (xo(p),yo{p)) into 
Eq. (jTT|> we obtain 



, , x 1, , pn/2 > Unrn 1 / -y H y (x ,yo) \ 



From this asymptotic expansion we obtain the logarith- mic correction term with coefficient of —1/2 exactly, for 
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any value of n. In fact, this asymptotic theory predicts 
that there exists such a logarithmic correction term with 
coefficient of —1/2 for a large class of lattice models when 
the two variables involved are proportional, that is, when 
the models are at fixed "density" . For those lattice mod- 
els which can be described by bivariate generating func- 
tions, this logarithmic correction term with coefficient of 
— 1/2 is universal when those models are at fixed "den- 
sity" . For the monomer-dimer model, this proportional 
relation is for s and to with s = pmn/2. An explicit 
calculation for n = 1 is shown in Appendix [5] 

For a fixed dimer density p and a fixed lattice width 
n, the first term of Eq. ^] is a constant and does not 
depend on to. we identify it as foo,n(p) 



(15) 



In all the following fitting experiments, we set £ 
-1/2 for Eq. 



IV. CYLINDER LATTICES 

For the monomer-dimer problem at a given dimer den- 
sity p in cylinder lattice strips, the sequence foo,n{p) con- 
verges very fast to f% (p) , as can be seen from a few sam- 
ple data in Table [I] In the table, values of /oo n(p) f° r 
p = 1/4, 1/2, 3/4, and 1 are listed. Two obvious features 
can be observed: (1) The function foa,n(p) is an increas- 
ing function of odd n, but a decreasing function of even 
n. Furthermore, for finite integer values of h and k, 



foo,2h{p) > h{p) > /oo,2fc+lG°)- 



(16) 



The value foo,n(p) oscillates around the limit value /2(p) 
from even n to odd n. (2) The smaller the value of p, the 
faster the rate of convergence of /oo,™(p) towards / 2 (/9)- 
Rational values of p are used for the calculations in Table 
[I] and no interpolation of a m ^ n (p) is used. The numbers 
of data points used in the fitting are listed in the paren- 
theses. 

As a check of the accuracy of the results, the data at 
p = 1 can be compared with the exact solution. For 
a cylinder lattice strip oo x n, the exact expression for 
/oo,«(l) reads as 



i/2 



/oo,n(l) = -InTT 



sin 



(2i- l)vr 



sin 



(2i- 1)tt 



(17) 

The results from Eq. El are listed as the last column 
in Table Q| As mentioned in the previous Section, all 
input data are exact integers and the logarithm of these 
integers is taken with high precision before the fitting. 
The only places where accuracy can be lost are in the 
fitting procedure as well as the approximation introduced 
by the fitting function Eq. [5] Comparisons of the data in 
the last two columns of Table [I] show that, as far as the 



fitting procedure is concerned, the calculation accuracy 
is up to 12 or 13 decimal place. 

Another check for the accuracy of the fitting proce- 
dure is through the exact expression Eq. IB4I of one- 
dimensional strip (n — 1) at various dimer density p. 
The data are listed in the first row of Table [I] By using 
Eq. I|15fl of the Pemantle and Wilson asymptotic method, 
we can also compare the fitting results with exact asymp- 
totic values for small values of n (data not shown). All 
these checks confirm consistently that the accuracy of 
the fitting procedure is up to 12 or 13 decimal place. See 
Section HXl for further discussions on this issue. 



The fast convergence of /oo,n(p) and the property of 
Eq. ^)]make it possible to obtain f2(p) quite accurately, 
especially when p is not too close to 1. Some of the values 
of /2(p) at rational p = p/q for small p and q are listed 
in Table ITT1 As in Table [I] no interpolation of a mj „(p) 
is used. The numbers in square brackets indicate the 
next digits for n = 16 (upper bound) and n = 17 (lower 
bound). The data show that when p < 0.65, the f2{p) 
is accurate up to at least 10 decimal place. It should 
be pointed out that the data listed are just raw data, 
showing digits that have already converged for n = 16 
and n = 17. If the pattern of convergence of these raw 
data is explored and extrapolation technique is used, as is 
done in Section lVll it is possible to get even more correct 
digits. As shown in Section IVTl the true value of fi{p) 
is not the average of fao,w{p) and foo,n{p)- Instead, it 
should lie closer to foo.nip)- 



LATTICES WITH FREE BOUNDARIES 



We also carry out similar calculations for lattice strips 
on to x n two-dimensional lattices with free boundaries, 
for n = 2, . . . , 16. A few sample data are shown in Table 
HTT1 In the table, values of f^ n (p) for p = 1/4, 1/2, 3/4, 
and 1 are listed. The data in Table II I II show that the 
sequence of f^. n (p) in lattices with free boundaries con- 
verges slower than that in cylinder lattices. Furthermore, 
in contrast to the situation in cylinder lattices, /^„(p) 
is an increasing function of n for < p < 1: f^ n {p) 
approaches f% (/?) (the same value as that for cylinder lat- 
•tices) monotonically from below. When p = 1, the func- 
tions /^ 2 fc(l) an d /<^2fc+i(l) are increasing functions, 
with /£ 2fe (l) > /S, 2fc _i(l) and /£ 2fe (l) > /£ 2&+1 (l). 
Due to the slow convergence rate and the lack of property 
like Eq. ^3 it is difficult to obtain / 2 (p) reliably from the 
data on the lattice strips with free boundaries. 

As we did in the previous Section, we also take advan- 
tage of the known exact solution for p — 1 as a check for 
the numerical accuracy of the fitting procedure. The ex- 
act result for lattice strips with free boundaries is given 



7 



TABLE I: The coefficient Co (/<x>,n(p)) for different n and p on cyiinder iattice strips oo x n. The numbers in parentheses are 
the number of data points used in the fitting. The first row for n = 1 is from the exact expression Eq. IB4I The last column is 
from the exact expression Eq. 1171 when p = 1. Rational p is used here and no interpolation of a m ,n(p) is used. 





1/4 


1/2 




3/4 




1 


1 


1 


0.358851778502358 


0.477385626221110 




0.4206:i2291*S07Sr, 




0.000000000000000 




1 


0.358851778501632 (113) 


0.477385626220963 


(226) 


0.420632291880650 


(113) 


3.6259082842339c-31 (451) 


0.000000000000000 


2 


0.443539035661245 (226) 


0.643863506776599 


(451) 


0.675072579831534 


(220) 


0.440686793509790 (901) 


0.440686793509772 


3 


0.441243226869578 (113) 


0.632058256526847 


(226) 


0.634554086596250 


(113) 


0.261133206162104 (451) 


0.261133206162069 


4 


0.441350608415009 (451) 


0.633331866235995 


(901) 


0.641840174628945 


(451) 


0.329239474231224 (901) 


0.329239474231204 


5 


0.441345086182334 (113) 


0.633177665529326 


(226) 


0.640045538037963 


(113) 


0.280932225367582 (451) 


0.280932225367553 


6 


0.441345392065621 (226) 


0.633198099780748 


(451) 


0.640485680552428 


(220) 


0.307299539523143 (901) 


0.307299539523125 


7 


0.441345374298049 (113) 


0.633195220523869 


(226) 


0.640363389854116 


(113) 


0.286180041989361 (451) 


0.286180041989328 


8 


0.441345375366775 (901) 


0.633195644943681 


(901) 


0.640398267527096 


(901) 


0.300105275372022 (901) 


0.300105275372003 


9 


0.441345375300735 (113) 


0.633195580174568 


(226) 


0.640387826199450 


(113) 


0.288315256713912 (451) 


0.288315256713877 


10 


0.441345375304906 (226) 


0.633195590329820 


(451) 


0.640391026472015 


(226) 


0.296935925720006 (901) 


0.296935925719986 


11 


0.441345375304640 (113) 


0.633195588702860 


(220) 


0.640390021971380 


(113) 


0.289391267149380 (451) 


0.289391267149350 


12 


0.441345375304658 (451) 


0.633195588968099 


(901) 


0.640390342494518 


(451) 


0.295260881552885 (901) 


0.295260881552868 


13 


0.441345375304658 (113) 


0.633195588924235 


(226) 


0.640390238745032 


(113) 


0.290008735546277 (451) 


0.290008735546247 


14 


0.441345375304656 (196) 


0.633195588931575 


(391) 


0.640390272712621 


(196) 


0.294265803657058 (781) 


0.294265803657028 


15 


0.441345375304652 (71) 


0.633195588930329 


(143) 


0.640390261482688 


(71) 


0.290395631458758 (285) 


0.290395631458698 


16 


0.441345375304640 (375) 


0.633195588930530 


(375) 


0.640390265226077 (375) 


0.293625491565320 (375) 


0.293625491565145 


17 


0.441345375304620 (28) 


0.633195588930470 


(57) 


0.640390263969286 


(28) 


0.290653983951606 (113) 


0.290653983951281 



TABLE II: List of /2(p) for different p. Numbers in square 
brackets indicate the next digits for n = 16 (upper bound) 
and n = 17 (lower bound). Rational p is used here and no 
interpolation of a mj „(p) is used. 



p 


Hp) 








1/20 


0.1334362263587 


1/10 


0.229899144084[8..9] 


3/20 


0.310823643168[1..2] 


1/5 


0.380638530252[1..2] 


1/4 


0.4413453753046 


3/10 


0.4940275921700 


1/3 


0.525010031447[5..6] 


7/20 


0.539305666744[5..6] 


2/5 


0.5775208675757 


9/20 


0.6088200746799 


1/2 


0.633195588930[4..5] 


11/20 


0.650499726669[5..8] 


3/5 


0.66044120984[2..5] 


13/20 


0.6625636470[2..4] 


2/3 


0.661425713[7..8] 


7/10 


0.65620036[0..1] 


3/4 


0.64039026[3..5] 


4/5 


0.6137181[3..4] 


17/20 


0.573983[2..3] 


9/10 


0.51739[1..2] 


19/20 


0.435[8..9] 


1 


0.29[0..3] 



VI. MAXIMUM OF FREE ENERGY AND THE 
MONOMER-DIMER CONSTANT 

It is well known that fd(p) is a continuous concave 
function of p and at certain dimer density p* , fd(p) 
reaches its maximum |37j . However, there is no ana- 
lytical knowledge of the location (p*) and value (fd(p*)) 
of the maximum for d > 1. As is shown in Appendix 
1X1 the maximum of fd(p) is equal to the monomer-dimer 
constant: fd(p*) = hd- Currently the best value for h 2 
is given in Ref. [H which gives h 2 = 0.6627989727 ± 
0.0000000001, with 9 correct digits. The location of the 
maximum, /?*, is controversial. Baxter gives the value of 
p* = 0.63812311 0, while Friedland and Peled state, 
"it is reasonable to assume that the value p* , for which 
X 2 (p*) = h 2 ,is fairly close to p(2) = ~ 0.6096118" 

(Here the original notation is used: p is our p and \ 2 {p) 
is our f 2 {p)) 

In this Section we use the same computational pro- 
cedure described in the previous sections to locate ac- 
curately the maximum of f 2 (p). Using rational dimer 
density p — p/q and choose appropriate p and q, we can 
locate the maximum to a fairly small region, as shown in 
Figure ^ 

With the interpolated data for a mi „(p), we can locate 
p* and f 2 (p* ) more accurately. As shown in Figures 
and |21 we find that 



by 



/£, n (l) = -ln 
n 



n 

i=l 



1 



1 



(18) 

The last column in Table lTTll lists the values given by Eq. 
ITH1 which can be compared with the calculated values 
from the fitting experiments in the column next to it. 
Again, as shown in the previous Section, the accuracy at 
p = 1 is up to 11 or 12 decimal place for most of the 
values of n. 



0.662798972831 < f 2 (p*) < 0.662798972845, 

where the value of /oo,n(p) for n — 16 is used as the upper 
bounds, and that for n — 17 as the lower bounds. From 
Figures [21 El and 0] we can locate p* as 

0.63812310 < p* < 0.63812312. 

The values of foo.n(p) around p* are listed in Table llVl 
Inspection of the convergent rate of these data for even 
and odd values of n suggests that for both sequences, the 
convergent rate is geometric. If we assume that 



foo.n(p) = h{p) ~ Urf 



(19) 
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TABLE III: The coefficient Co (/oo,n(p)) f° r different n and p on lattice strips oo x n with free boundaries. The numbers in 
parentheses are the number of data points used in the fitting. The last column is from the exact expression Eq. I18l when p = 1. 
Rational p is used here and no interpolation of a m ,n(p) is used. 



1/4 



1/2 



3/4 



1 



1 



0.406768721898144 
0.418805029581931 
0.424677898377694 
0.428121453697918 
0.430386238446347 
0.431988836086132 
0.433182588323077 
0.434106231271596 
0.434842114982272 
0.435442204742419 
0.435940910510948 
0.436361922501114 
0.436722083762241 
0.437033697160078 
0.437305958542365 



(101) 

(50) 

(201) 

(50) 

(101) 

(50) 

(401) 

(50) 

(101) 

(50) 

(201) 

(39) 

(26) 

(13) 

(44) 



0.567460205873414 
0.589202338098224 
0.600481083876114 
0.607125402184205 
0.611530695170404 
0.614662382427737 
0.617003233867274 
0.618819152717284 
0.620268892851619 
0.621453058127923 
0.622438494443121 
0.623271352033866 
0.623984518924941 
0.624602065315795 
0.625142013068189 



(201) 

(101) 

(401) 

(101) 

(201) 

(101) 

(401) 

(101) 

(201) 

(101) 

(401) 

(78) 

(51) 

(26) 

(44) 



0.55061882 12751,01) 
0.577814537070212 
0.593860234282314 
0.603150396985283 
0.609386552832152 
0.613828224552787 
0.617157805951205 
0.619745522952440 
0.621814606701445 
0.623506740930563 
0.624916337018327 
0.626108703477498 
0.627130461956123 
0.628015783739589 
0.628790285827699 



(101) 

(50) 

(201) 

(50) 

(101) 

(50) 

(401) 

(50) 

(101) 

(50) 

(201) 

(39) 

(26) 

(13) 

(44) 



0.240605912529824 
0.219492982820793 
0.260998208772619 
0.252922288709197 
0.269862305348313 
0.265557149993036 
0.274751610011806 
0.272072662436541 
0.277844105572086 
0.276016066623932 
0.279975752031904 
0.278648778924217 
0.281534000787684 
0.280526932170974 
0.282722754819597 



(401) 

(201) 

(401) 

(201) 

(401) 

(201) 

(401) 

(201) 

(401) 

(201) 

(401) 

(155) 

(101) 

(51) 

(44) 



0.240605912520802 
0.219492982820803 
0.260998208772539 
0.252922288709162 
0.269862305348238 
0.265557149992917 
0.274751610011700 
0.272072662436436 
0.277844105571997 
0.276016066623911 
0.279975752031819 
0.278648778924210 
0.281534000780413 
0.280526932164772 
0.282722752409010 



-S 0.66278 




0.66279897290 



0.66279897285 



0.66279897280 



0.66279897275 



0.66279897270 



0.66279897265 



0.66279897260 



0.66279897255 



0.66279897250 




FIG. 1: (Color online) The function of f2(p) in the region 
of 19/20 < p < 9/14. All data points use rational p, so no 
interpolation is used. 



FIG. 2: (Color online) The function of foo,n(p) in the region 
of 0.638110 < p < 0.638130 for n = 16 (upper curve) and 
n — 17 (lower curve). 



then the data points at p = 0.63812311 of n = 12, 
14, and 16 can be used to obtain an extrapolated value 
of f 2 (p) = 0.6627989728336, while the data points of 
n = 13, 15, and 17 give another extrapolation value 
f 2 (p) = 0.6627989728341. Together these two extrapo- 
lation values converge to f 2 (p*) = 0.662798972834, with 
11 correct digits. 

We can also get the same conclusion graphically from 
Figure |3J By inspecting the pattern of the data points 
of different values of n in the inset of Figure |31 we notice 
that the difference between the data points of n = 14 and 
ri = 16 is bigger than the difference between n = 15 and 
n = 17. This indicates that the true value of f 2 {p*) ncs 
closer to the data point of n = 17 than the data point of 
n = 16. From Figure [3] we are quite sure that the 11-th 
digit of f 2 (p*) is 3 instead of 4, and the 12-th digit is 
probably 4, as indicated by the two extrapolation values 
mentioned above. 

The value of . hip* ) is in excellent agreement with that 
reported in Ref. [l^, which gives 9 correct digits (Fried- 
land and Peled also guess correctly the 10-th digits as 8). 
The value also agre es with that in Ref. H^. which gives 
8 correct digits 12]. The value of p* is exactly that of 



Baxter [Fjj. 

By using field theoretical method, Samuel uses the fol- 
lowing relation to transform the activity x into a new 
variable to [lfj| 



(1 -Auf 



(20) 



This relation is very close to the one used by Nagle 
j8|- By substituting this relation into Gaunt 's series 
expansions Samuel obtained new series for various 
lattices, including the rectangular lattice (Eq. (5.12) 
of Ref. |I3). The value of monomer-dimer constant in 
two-dimensional rectangular lattice can be calculated at 
x = 1 by using his series as 0.66279914. As we can see, 
this only gives five correct digits. Nagle used the follow- 
ing transform § 



(21) 



(l-3w) 2 



By using Gaunt's series, a value of 0.6627988 is obtained, 
with six correct digits. 

To conclude this section, we compare our results on 
the maximum of f 2 (p) with the approximate formulas of 
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TABLE IV: The coefficient Co (/oo,n(p)) for different n and p on cylinder lattice strips oo x n around p* . The numbers in 
parentheses are the number of data points used in the fitting. 



0.63812309 0.63812310 0.63812311 0.63812312 



1 





470643631091106 


(880) 





470643628559868 


(880) 





470643626028629 


(880) 





470643623497390 


(880) 


2 





683451694063943 


(901) 





683451695019491 


(901) 





683451695975038 


(901) 





683451696930585 


(901) 


3 





659839104062378 


(901) 





659839103873019 


(901) 





659839103683659 


(901) 





659839103494298 


(901) 


4 





663319985040839 


(901) 





663319985089007 


(901) 





663319985137175 


(901) 





663319985185343 


(901) 


5 





662701144811933 


(901) 





662701144800592 


(901) 





662701144789251 


(901) 





662701144777910 


(901) 


6 





662818978977777 


(901) 





662818978980627 


(901) 





662818978983477 


(901) 





662818978986327 


(901) 


7 





662794695257766 


(901) 





662794695257048 


(901) 





662794695256327 


(901) 





662794695255611 


(901) 


8 





662799924786436 


(901) 





662799924786622 


(901) 





662799924786807 


(901) 





662799924786992 


(901) 


9 





662798754939549 


(901) 





662798754939501 


(901) 





662798754939453 


(901) 





662798754939404 


(901) 


10 





662799023857733 


(901) 





662799023857746 


(901) 





662799023857758 


(901) 





662799023857771 


(901) 


11 





662798960670265 


(901) 





662798960670262 


(901) 





662798960670259 


(901) 





662798960670256 


(901) 


12 





662798975775941 


(901) 





662798975775943 


(901) 





662798975775944 


(901) 





662798975775944 


(901) 


13 





662798972113454 


(901) 





662798972113453 


(901) 





662798972113445 


(901) 





662798972113451 


(901) 


14 





662798973011855 


(781) 





662798973011855 


(781) 





662798973011855 


(781) 





662798973011855 


(781) 


15 





662798972789303 


(570) 





662798972789304 


(570) 





662798972789304 


(570) 





662798972789303 


(570) 


16 





662798972844882 


(375) 





662798972844882 


(375) 





662798972844883 


(375) 





662798972844883 


(375) 


17 





662798972830869 


(226) 





662798972830871 


(226) 





662798972830870 


(226) 





662798972830869 


(226) 



0.662798972850 



0.662798972845 ■ 



0.662798972840 



0.662798972835 - 



n=16 — |— 

n=i7 — X-- 



0.66279897305 



0.66279897283 



□□□□□□□□□□□□□□□□□□□ 

0.66279897275 

0.638122 0.638123 0.638124 

0.662798972830 - v X -. X " X " * ' * " ^ -X-X-X-- X - V 



0.662798972825 



Chang [ljj and Lin and Lai |l8j . The Chang's approxi- 
mate formula is given below 



/ 2 c (p) 



i 



[4 In 4- (4-p) ]n{4-p)+p In p+2{l-p) ln(l-p)-2/5 In 4] , 



which gives p* = 0.634641 and f 2 {p*) = 0.661355. The 
approximate formula of Lin and Lai is 

J 2 LL (P) = -| m | - 0.9030(1 - p) ln(l - p) - 0.05645p 

which gives p* = 0.638057 and f 2 (p*) = 0.66282235. Al- 
though the two formulas are quite simple, they give ef- 
fective approximation with respect to p* and f 2 (p*)- 



FIG. 3: (Color online) The function of /oo,n(p) in the region 
of 0.6381221 < p < 0.6381239 for n = 16 (upper curve) and 
n = 17 (lower curve). In the inset the data points from n = 14 
and n = 15 are also shown. From top to bottom in the inset: 
n = 14, n = 16, n = 17, and n = 15. 



0.662798972844884 
0.662798972844882 
0.662798972844880 
0.662798972844878 
0.662798972844876 
0.662798972844874 
0.662798972844872 
0.662798972844870 
0.662798972844868 
0.662798972844866 
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+ 
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+ 




+ 

+ - 


+ 







FIG. 4: (Color online) The function of foo,n(p) in the region 
of 0.63812301 < p < 0.63812319 for n = 16. 



VII. COMPARISON WITH BAXTER'S 
RESULTS 

Using variational approach, Baxter calculated h 2 (x) 
(which is In k using his notation) and 9(x) (which is 2p 
by his notation) for different values of dimcr activity x 
(s 2 by his notation). By using Eqs. IA4I and IA5I we can 
compare our results with Baxter's results in his Table II. 
For each of his data point at a dimer activity x, we cal- 
culate f2(p) with p = 8(x). Then his h 2 (x) is converted 
to f 2 (p) = h 2 (x) — | ln(x). The comparisons are shown 
in Table It should be pointed out that in Baxter's 
data, extrapolation is used for the sequence to obtain 
9(x) and h 2 {x) when x~ 1 / 2 is small (x^ 1 / 2 < 0.3 for 
h 2 (x) and x~ 1 ' 2 < 0.5 for 0(x)), while no extrapolation 
is used in our data: we only look at the digits that have 
been converged for n = 16 and n = 17. Although the ex- 
trapolation used in Baxter's data makes the comparison 
less direct, we still see that the agreement is excellent. It 
seems that Baxter's method converges faster for p very 
close to 1 (again the extrapolation factor has to be con- 
sidered) , and our method is more accurate when p is not 
too close to 1. As in Section llVl we only present the raw 
data here. If extrapolation is used, more correct digits 
can be obtained. 
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TABLE V: Comparison with Baxter's results. Numbers in 
square brackets indicate the next digits for n = 16 (upper 
bound) and n = 17 (lower bound). 



-1/2 



h0>) 



fi{p) 



0.00 
0.02 
0.05 
0.10 
0.20 
0.30 
0.40 
0.50 
0.60 
0.80 
1.00 
1.50 
2.00 
2.50 
3.00 
3.50 
4.00 
4.50 
5.00 



1.0 

0.994176 

0.9836216 

0.96456376 

0.924706050 

0.8846581140 

0.8453815864 

0.8072764728 

0.7705280966 

0.7013863228 

0.6381231092 

0.5042633294 

0.4006451804 

0.3211782498 

0.2603068980 

0.2134739142 

0.17715243204 

0.14869898092 

0.126162903820 



0.29[0..3] 

0.319[2..8] 

0.355[0..2] 

0.4047[5..8] 

0.4810[8..9] 

0.536892(1. .4] 

0.5782845[2..9] 

0.608814[3..4] 

0.63085609[6..8] 

0.655894637[3..5] 

0.6627989728[3..4] 

0.6349499289380[4..9] 

0.5779686472227(1. .4] 

0.5140847735884[4..6] 

0.4528361791290[7..9] 

0.3978378948658[1..3] 

0.3499573614350[0..2] 

0.3088705309099[2..6] 

0.273811439807[1..2] 



0.291557 

0.3194631 

0.35510683 

0.404771005 

0.4810887477 

0.5368922350 

0.5782845477 

0.6088143934 

0.6308560970 

0.6558946374 

0.6627989726 

0.6349499290 

0.5779686472 

0.5140847737 

0.4528361790 

0.3978378949 

0.3499573615 

0.3088705306 

0.2738114398 



VIII. HIGH DIMER DENSITY NEAR CLOSE 
PACKING 

It is well known that phase transition for the monomer- 
dimer model only occurs at p = 1 0]. Since the close- 
packed dimer system is at the critical point, it is interest- 
ing to investigate the behavior of the model when p — > 1. 
Using the similar computational procedure outlined be- 
fore, the following results are obtained at high dimer den- 
sity limit: 



/oo,n(p) ~ /™(!) + 



-(1 — p) ln(l — p) n is odd 
, -5(1 — p) — p) n is even 

(22) 

where /<^ t J l lcc (l) is the free energy of close-packed lat- 
tice with width n, and is given, based on the boundary 
condition, by Eq. El (cylinder lattices) or Eq. ^| (lat- 
tices with free boundary condition). Eq. 1221 for n = 1 is 
verified from the exact result as shown in Eq. IB5I The 
result is also confirmed for other values of n by using the 
Pemantle and Wilson asymptotic methods for multivari- 
ate generating function, as described in Section ITTT1 For 
space limitation these confirmations are not presented in 
this paper. 

The dependence of the asymptotic form of /oo,n(p) on 
the parity of the lattice width n as shown in Eq. [22 re- 
minds us of the results reported previously for monomer- 
dimer model with fixed number of monomers v [29|, in 
which the coefficient of the logarithmic correction term of 
the free energy depends on the parity of the lattice width 
n. These two results are consistent with each other. If 
we substitute the relation v = (1 — p)mn into Eq. 1221 we 
will get the logarithmic correction term with coefficient 




0.999 0.9991 0.9992 0.9993 0.9994 0.9995 0.9996 0.9997 0.9 
P 

FIG. 5: (Color online) Fitting result for /2(p) as p — > 1. 



v for odd n, and v/2 for even n. More discussions about 
this asymptotic form will be found in Section IIXI (Eq. 

mi. 

We also investigate the behavior of f%{p) (for infinite 
lattice) as p — *■ 1. Since foo,n(p) does not converge fast 
enough as p — > 1 (Table HJ, we use weighted average of 
foo,ie(p) an d foo,n{p) as an approximation of /2(p)- The 
weights are calculated from the exact results at p = 1. 
Fitting these data to the following function 

h{p) = G/tt + 2(1 - p) In (1 - p) + h(l - p), (23) 

we obtain 7 w 1.69775 and bi w 0.427832. No other 
reasonable form of functions other than Eq. I|23|l gives 
better fit. Including a term of (1 — p) 2 in Eq. 1231 leads 
to only slight changes in the values of 7 and b\ The data 
and the fitting result are shown in Figure 

Using the equivalence between statistical ensembles 
discussed in Appendix^ we can relate our results with 
Gaunt's series expansions |(|- Plugging /2G0) as in Eq. 
OS into f(p) + § ln(x) (see Eq. |MJj, differentiating with 
respect to p, and solving for p, we obtain the average 
dimer density 9{x) at the activity x. Expressing x as a 
function of 8, we have 



(1 - eyi 



This is in the same form of Eq. (3.7) of Gaunt [fj. If 
we put in the values of 7 and b\, we can estimate the 
amplitude A = exp(26i — 7) = 0.4308. Gaunt obtains 
through series expansions 7 = 1.73 ±4 and A = 0.3030 ± 
4, and conjectures that 7 = 7/4. Our results support the 
same functional form, and the numerical values are close 
to these obtained by Gaunt's series analysis. As for the 
conjectured value of 7, the current data seem to indicate 
a value lower than 7/4. In fact, the data presented here as 
well as theoretical arguments (not shown here) indicate 
that 7 = 5/3. More discussion on this constant can be 
found in the next Section. 
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IX. DISCUSSION 

In Section [illl we show by computational methods that 
there is a logarithmic correction term in the free energy 
with a coefficient of —1/2. By introducing the newly 
developed Pemantle and Wilson asymptotic method, 
we give a theoretical explanation of this term. We 
also demonstrate that this term is not unique to the 
monomer-dimer model. Many statistical lattice models 
can be casted in the form of bivariate (or multivariate) 
generating functions, and when the two variables are pro- 
portional to each other so that the system is at a fixed 
"density" , we will expect such a universal logarithmic 
correction term with coefficient of —1/2. We anticipate 
more applications of this asymptotic method in statisti- 
cal physics in the future. 

The Pemantle and Wilson asymptotic method not only 
gives a nice explanation of the logarithmic correction 
term and its coefficient found by computational means, 
but also gives exact numerical values of foo,n(p) for small 
n (the width of the lattice strips). These exact values 
can be used to check the accuracy of the computational 
method, as already mentioned in Section llVI In Section 
IIIII we discuss how this can be done. The denominator 
H(x, y) of the bivariate generating functions is derived 
from the characteristic function of the matrix M n , and 
the size of M n is given by u c [n) in Section[n]for cylinder 
lattices, and in Ref. |2^ for lattices with free boundaries. 
For small n, the size of M n is small enough so that the 
characteristic function can be calculated symbolically. As 
n increases, however, the size of M n increases exponen- 
tially: u c (l7) = 4112 for cylinder lattice when n = 17, 
and itfb(16) = 32896 when n = 16. It is currently imprac- 
tical to calculate the characteristic functions symbolically 
from matrices of such sizes to get H(x, y) of the corre- 



sponding bivariate generating functions, so the Pemantle 
and Wilson method cannot be applied when the width of 
the lattice n becomes larger. Even when H(x, y) is avail- 
able, it is of the order of thousands or higher, which will 
lead to instabilities in the numerical calculations. The 
computational method utilized here, however, can still 
give important and accurate data in these situations. 

Previously we demonstrated that when the monomer 
number v or the dimer number s are fixed, there is also 
a logarithmic correction term in the free energy |28l l29j | . 
When the number of dimers is fixed (low dimer density 
limit), the coefficient of this term equal to the number 
of dimers. When the number of monomers is fixed (high 
dimer density limit), the coefficient, however, depends on 
the parity of the lattice width n: it equals to v when n is 
odd, and v/2 when n is even. In this high dimer density 
limit, as m — > oo, dimer density p — > 1. In this paper we 
focus on the situation where the dimer density is fixed, 
and find that again there is a logarithmic correction term, 
but this time its coefficient equals to —1/2 and does not 
depend on the parity of the lattice width. Why does the 
dependence of the coefficient on the parity of the lattice 
width disappear as dimer density p — > 1 and the lattice 
becomes almost completely covered by the dimers? 

This seemingly paradoxical phenomenon can be ex- 
plained as follows. When the number of monomer v is 
fixed and as m — > oo, if we can put v = (1 — p)mn into 
Eq. then the term of (1—p) ln(l— p) leads to a term of 
win to/ '(mn) when n is odd, and a term of vlnm/(2mn) 
when n is even. At the same time, the logarithmic cor- 
rection term with —1/2 as coefficient (— lnm/(2mn), the 
second term in Eq. 1141) , gets canceled out by the a term 
of — ln(l — p)/(2mn) from the third term in Eq. 1141 as 
to — ► oo and p — > 1. Putting Eqs. 1141 and l2~2l together, 
we have for finite n, as m — > oo and p — ► 1, 



J 



fm,n(p) 



/^ C0 (l) 



-(l-p)ln(l-p) 
-I(l-p)ln(l-p) 



2TO7J 



■ In i 



2mn 



ln(l - p) 



n is odd 



n is even 



(24) 



r 



This expression can be checked with the explicit formu- 
las for one-dimensional lattice, Eqs. IB2I and IB5I By 
using the relation v = (1 — p)mn, we see from Eq. 1241 
that as to — > oo, when the monomer number v is fixed, 
the dependence of the logarithmic correction term on the 
parity of n comes from the second term in the equation; 
the third and fourth terms cancel each other out. On the 
other hand when the dimer density p is fixed, the only 
logarithmic correction term comes from the third term of 
Eq. EH with coefficient -1/2. 

As n — ► oo, f%(p) also has a term of (1 — p) ln(l — p) 
(Eq. \2'6\ . whose coefficient is estimated as —0.85 (Section 
|VII1*)| . This value is closer to -5/6 = -0.83 than to the 
conjectured value —7/8 = —0.875 by Gaunt Runnels' 



result, however, seems to be closer to Gaunt's result |32| . 
It should be recognized that, as pointed out previously 
H> M El H3 as well as in the present work, the con- 
vergence is poorest when the dimer density is near close 
packing. On the other hand, theoretical calculations un- 
derway (not shown here due to space limitation) indeed 
indicate that this coefficient of (1 — p) ln(l — p) for the 
infinite lattice is —5/6, or equivalently, 7 = 5/3. 

It is well known that there is a one to one correspon- 
dence between the Ising model in a rectangular lattice 
with zero magnetic field and a fully packed dimer model 
in a decorated lattice [23, ■ By using a similar method 
|l4j , it has been shown that the Ising model in a non-zero 
magnetic field, a well-known unsolved problem in sta- 
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tistical mechanics, can be mapped to a monomer-dimer 
model with dimer density p < 1. The investigation of 
the monomer-dimer model near close packing is of inter- 
est within this context. 

As mentioned in Section IVTl several authors have ap- 
plied field theoretic methods to analyze the monomer- 
dimer model (for example, Ref. ITol) . In such studies, the 
monomer-dimer problem is expressed as a fermionic field 
theory. For close-packed dimer model, the expression is 
a free field theory with quadratic action, which is exactly 
solvable as expected. For general monomer-dimer model, 
the expression is an interacting field theory with a quartic 
interaction term, and self-consistent Hartree approxima- 
tion is used to improve the Feynman rules to derive the 
series expansions. The transforms that are obtained us- 
ing these methods (such as Eq. [2U1 which is similar to 
Eq. I21f> make the series expansions converge in the full 
range of the dimer activity. The accuracy of these cal- 
culations, however, is not comparable to the accuracy of 
the computational method reported here, possibly due to 
the limited length of the series expansion. 



APPENDIX A: EQUIVALENCE OF 
STATISTICAL ENSEMBLES 



The connection for the average dimer coverage can also 
be obtained by using Eqs. 12*1 and |A"T1 as 



6(x) = lim m , n (x) = lim — d ^ n ^ mn ^ = p *( x ) 

m,n^oo m, oc mn O m X 

(A4) 

with the understanding that at p* , F 2 (p, x) = f 2 (p) + 
f \nx, not f 2 (p), reaches its maximum. Substituting Eq. 
4] into Eq. EH we obtain 



h 2 (x) = f 2 (9(x)) + 



9{x) 



In. 



(A5) 



In Section IVII the excellent agreement of our result 
of h(p*) with the result on h 2 of Friedland and Peled 
[l2( has already been demonstrated. In Ref. H2l what 
is calculated is actually h 2 - Eq. IA3I makes it possible to 
compare our result with that of Friedland and Peled. Eq. 
IA3l is proved as a theorem for the specific monomer-dimer 
problem in Ref. ^3 as Theorem 4.1. 

Since there is an analytical solution to the one- 
dimensional lattice problem, Eqs. IA1I and IA4I can be 
confirmed for the one-dimensional lattice by explicit cal- 
culations, as shown in Appendix iBl 



Throughout the paper our focus is on the functions 
fm,n(p), foo,n(p), or h{p) at a given dimer density p. 
These functions are in essence properties of the canoni- 
cal ensemble. In this Appendix we make the connection 
between f 2 (p) and the functions of 8{x) and li2(x) as de- 
fined in Eqs. 0and0] which are properties of the grand 
canonical ensemble. The results are used in Section lYlHI 
to compare the results of near close packing dimer density 
with Gaunt's series analysis, and in Section IVIII to com- 
pare our results with those of Baxter, whose calculations 
are carried out in terms of 9{x) and h 2 (x) |l9j . 

Suppose at p — p* , the summand a m ^ n (p)x mnp ^ 2 in 
Eq. El reaches its maximum. By using the standard ther- 
modynamic equivalence between different statistical en- 
sembles [13, Chap. 4], we have 



h,2{x) = lim 



h\Z m ^ n (x) 



m,n^oc mn 



= lim 



m,n — >oo 



lim 

m,n — >oo 



ln E <p<l a r»,n(p)z 

mn 

lna m , n {p*)x mn P*/ 2 



mnp/2 



Dili 



= h(p*) + y lnx. 
In other words, if we define F 2 (p, x) = f 2 (p)+o l na; ; then 



(Al) 



APPENDIX B: EXPLICIT RESULTS FOR 
ONE-DIMENSIONAL LATTICE 

In this Appendix we summarize some exact results for 
the one-dimensional lattice (n = 1) which are useful to 
compare and check the results for lattices with width 
n > 1. When n = 1, the problem is a special case of the 
so-called "parking problem" in one-dimensional lattice 
in which a fc-mer covers k consecutive lattice sites in a 
non-overlapping way. Various methods exist which lead 
to closed form solutions to the general case of interacting 
fc-mers (for example, see Refs. 133. l4(i[) . For the monomer- 
dimer model, k = 2 and there is no interaction between 
the dimcrs. The number of ways to put s dimers in the 
m x 1 lattice is known as 



m — s 
s 



(Bl) 



From this expression, in the next subsection we derive the 
asymptotic expression of the free energy by using the tra- 
ditional method. As an illustration, later we also give an 
explicit demonstration of Pemantle and Wilson's asymp- 
totic method as it is applied to the bivariate generating 
function. 



h 2 {x) = max (f 2 {p) + - In a:) = max F 2 (p, x). (A2) 

0<p<l Z 0<p<l 

As a special case, the monomer-dimer constant is the 
maximum of the function f 2 (p) by setting x = 1 

h 2 = max / 2 (p) = / 2 (p*). (A3) 

0<p<l 



Canonical ensemble 



From the explicit expression of Eq. IB1I we can get the 
asymptotic expression of the free energy by using Stirling 
formula when < p < 1: 
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lna m ,i(s) , 1 1 2 - p 
fm,l{P) = = foo,l{P) ~ TT~ m TO + 7T- m "Ti V 

to 2m 2m p(l — p) 



1 



1 



(2-P) 



2j-l 



^i(2j-l)m2j 



1 



p 2i-i 2 2 J- 1 (1-/j) 



2i-i 



E 



i 



2 2i-2 B , 



2.7 



j(2j-l) V(2-P) 



1 



1 



1 



pjy-i 2 2 J-!(1 -p) 2 ^ 1 
I 



2(2j - 1) 



E 



4jm 



2j+l 



(B2) 



(B3) 



where 

/=o,i(p) = (l-f)Ml-f)-f lnf-(l-p)ln(l-p) (B4) 

and i?2j are the Bernoulli numbers. 

From Eqs. IB2l or lB3l it is evident that for n = 1, the 
coefficient of the logarithmic correction term is I = — 1/2 
for < p < 1. 

The asymptotic expression of /oo,i(p) (Eq. IB4(1 at p = 
1 is given by 



To calculate 9(x) using Eq. we need to evaluate the 
sum 

m/2 



S(m) = 



fooAp) ~ -(1-p) m(l-p)-(ln2-l)(l-p)- 



(1-P) 



2i+l 



^ 2i(2z+l) 
(B5) 

From this asymptotic expression we can see that the co- 
efficient of (1 — p) ln(l — p) is exactly —1, as in Eq. 1221 for 
odd values of n. By combining Eqs. IB 21 and IB 51 together 
we confirm Eq. 1241 for n = 1 at high dimer density limit. 



2. Grand canonical ensemble 

In this section we calculate various quantities associ- 
ated with the grand canonical ensemble. The configura- 
tional grand canonical partition function (Eq. ^) is 



m/2 

■y - / III 



To get a closed form of Z mi i(x), we use the WZ method 
(Wilf-Zeilberger) ^l| to get the following recurrence of 

Z m ,i( x ) 

xZ m% \{x) + Z m+1 .i(x) - Z m+2t i(x) = 0, 
from which we obtain the closed form solution as 

Z m>1 (x) = -j=L= - (3™ +1 ] (B6) 

VI + 4a; 

where 



1± V1 + 4.T 



s=0 

Again by using WZ method, we obtain the following re- 
currence for S(m) 

(m + 2)xS(m) + (m + l)S(m + 1) - mS(m + 2) = 0. 

To solve this recurrence, we use the generating function 
of S(m): Gs{z) = J2 m S{m)z m , and get Gs{z) from the 
recurrence as 



Gs(z) 



(1 



xz^ 



From Gs{z) a closed form expression of S(m) can be 
found as 



5(m) 



1 + 4x 



1 



=)/?r + (™+ 



i 



VT + 4x 1 VI +4x 

Substituting Eqs. IB6I and lB7l into Eq. we obtain 

1 



1 



y/T+Tx 

Using Eq. IB6l we can calculate h\(x) as 



(B7) 
(B8) 



fc l( x)= lim ln ^ l(x) =lni±4±^ (B9) 



It is known that there are multiple methods to solve the 
one-dimensional lattice model. For example, transform 
matrix method can also be used |4^. In this case, the 
transform matrix is 



x 

1 1 



whose eigenvalues are Ai,2 = (1 ± \/l+Tx)/2. As to — > 
°° 5 Z m ,i{ x ) ~ so we obtain /ii(x) as Eq. IB 91 From 
0i (x) = 2<91n Ai/<91nx we obtain Eq. IB8I 
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3. Equivalence of statistical ensembles 

The confirmation of the equivalence of ensembles for 
the special case of x = 1 (Eq. IA3|) has been done in Ref. 
IT2I Here we carry out the explicit calculations for the 
general case of arbitrary dimer activity x. 

If we take the derivative of function Fi (p, x) = fi (p) + 
fin a;, where fi(p) = foo,i{p) is given in Eq. IB4I solve 
for p, and retain only the solution in [0, 1], we have 



1 



1 



s/T+4x 



which is the same as Eq. IB8I If we put the value of p* 
into Fi (p, x) , we obtain the maximum of F\ (p* , x) : 



mation of Eq. IB6I or by using the characteristic function 
of Mi, or by using Eq. (23) of Ref. H] (by setting the 
interaction parameter a = 1 and the size of fc-mer as 2), 
to get 



Fi(x,y) = 



1 - v - xy 



(BIO) 



Here G(x,y) — 1 and H(x,y) = 1 — y — xy 2 . Solving 
the two equations in Eq. ^| we get (xo,yo) as (xq — 
s(m — s)/(m~ 2s) 2 , yo — (m — 2s)/(m— s)). Substituting 
(xq, yo) hito Eq. ^2 we obtain 



, f ( , p. , , 1 + V1 + 4.T 
max (fx (p) + -lnx)=]n 

0<p<l i z 



= hx(x) 



j „ ^(m~ S ) m - s+1 / 2 (m-2 S )- m+2s - 1 / 2 ( S )- s - 1 / 2 . 
v 2n 



4. Application of Pemantle and Wilson asymptotic 
method to the bivariate generating function 

The bivariate generating function of Eq, 1131 can also 
be obtained in multiple ways, for example by direct sum- 



By putting s = pm/2, the first three terms of Eq. IB2l are 
recovered, including the term of logarithmic correction 
— him / (2m) . Higher order terms can also be obtained if 
more terms of the asymptotic expressions are used 26] . 
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